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1.   Introduction 

Consider  the  following  inverse  eigenvalue  problem: 

Given  n  real  symmetric  n:<n  matrices  A,  ,  .  .  .  ,A   and  n  eigen- 
*         *  "  . 

values  X,  <...  <X  ,  find  the  coefficients  c, , . . .c      such  that  the  matrix 
1  —    —  n  In 

A(c)  =  c^A^+. . .+c^A^  (1.1) 

T 
has  the  given  eigenvalues.   Here   c  =  Ic,,...,c  ] 

A  slight  reformulation  of  (1.1)  describes  the  additive  inverse 
eigenvalue  problem,  which  arises  in  the  solution  of  inverse  Sturm- 
Liouville  problems.   In  practice  it  happens  frequently  that  only  some 
eigenvalues  are  given.   However,  for  the  purpose  of  analysis  it  is 
convenient  to  consider  the  case  where  the  number  of  parameters,  eigen- 
values and  the  order  of  the  matrices  is  the  same. 

Another  area  where  problem  (1.1)  arises  is  in  shell  model  com- 
putations in  nuclear  spectroscopy  (see  Brussard  and  Glaudemans  (1977)). 
There  A  is  the  Hamiltonian,  the  variables  {c^}  are  the  interactions  of 
one  and  two  bodies,  and  the  matrices  {A.}  represent  the  result  of 
adding  and  s^TTunetrizing  (or  antisymmetrizing)  the  effect  of  many 
particles. 

We  will  now  briefly  describe  an  inverse  Sturm-Liouville  problem. 
Consider  the  boundary  value  problem 

-  u"  +  p(x)u  =  ,\u      on  [0,7t]  (1.2) 

with  some  boundary  conditions,  for  example   u(0)  =  niir)    =    0.   Suppose 

*    DO 

that  p(x)  is  unkno'A^.   Given  a  spectrum.  {a^;-|_  ^or    the  probler,  (1.2), 
can  we  determine  p(x)?   Let  us  discretize  (1.2):   let  h  =  ^^  ,  u^  =  u  (ih)  , 
p.  =  p(ih)  ,  i  =  1,  .  .  .  ,n. 
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Then,  using  finite  differences  to  approximate  u", 


u. 


'i-1 


u 


0  . 


•O     n+1 
In  matrix  form  (1.3)  can  be  written  as 


(i=l,...,n) 


(1.3] 


where 


^0'? 


2-1 

-1  2  -1 
-1  . 


D  =  diag  (p^ ) 


(1.4) 


The  problem  is  therefore:   given  A-  and  a  spectrum  {Aj^}-.  find  a  diagonal 
matrix  D  such  that  A_+D  has  the  given  spectrum.   This  is  the  additive 
inverse  eigenvalue  problem.   It  can  be  written  in  the  form 


A(c) 


y    c.A. 
i4i    ^  1 


(1.5) 


with  A.  =  e.e.,  where  e.  is  the  i-th  column  of  the  identity  matrix. 
1111 

The  numerical  methods  and  the  convergence  results  that  we  give  for  the 

form  (1.1)  also  apply  to  (1.5). 

The  question  that  arises  immediately  is:   when  does  there  exist 

a  solution  to  problem  (1.1)  or  (1.2),  and  when  is  it  unique?   The  answer 

has  been  given  for  various  special  cases.   We  refer  the  reader  to  Borg 

(1946),  Hadeler  (1968),  Hald  (1972)  and  Hochstadt  (1974)  for  some  basic 

results  and  more  references.   In  the  case  of  the  discrete  prob.lem  (1.3) 


one  can  show  that  if  p  is  symmetric  (p. 


.  n 


^n+1 


^)    then  the  set  {X^]^ 


determines  {?■}"-,  •   If  p  is  not  symmetric  one  must  provide  two  sets  of 
eigenvalues  corresponding  to  two  sets  of  boundary  conditions. 

In  this  paper  we  will  assume  that  (1.1)  has  a  solution  and  will 
concentrate  on  how  to  compute  it  numerically. 


2.   Formulation  of  the  Numerical  Methods 

V7e  will  consider  only  the  problem  (1.1).   We  note  that  stable 
algorithms  exist  for  some  special  cases  of  (1.1);  see  de  Boor  and 
Golub  (1978)  ,  which  treats  the  construction  of  a  Jacobi  matrix  from 
"ictral  data. 


Let  '•  ,  (c)  <  .  .  .  <  X  (c)  be  the  eigenvalues  of  A  (c)  =   /_   c  •  A.  .   W^ 
1     —    —  n  ^_-^       1  1 

will  consider  two  formulations  of  the  problem. 


Forr.ulation  I.   Solve  the  nonlinear  system. 


f  (c] 


A  (c)  -  X 


=  0 


(2.5] 


Formulation  II.   Solve  the  nonlinear  system 


g(c)  = 


det(A(c)  -  X^I) 


det(A(c)  -  X^I) 


(2.6) 


The  first  formulation  is" the  most  natural  one  and  has  been  used  by 
Downing  and  Householder  U956)  in  the  case  of  the  additive  inverse 
eigenvalue  problem,  and  by  Kublanovska ja  (1970)  for  the  general  linear 
form  (1.1)  The  second  formulation  has  been  proposed  by  Biegler-Konig 
(1981) .   We  will  give  some  arguments  that  indicate  that  the  latter  has 
some  major  disadvantages. 

Consider  for  example  the  2x2  case  where  the  {A.}  are  diagonal. 
Then 


A(c) 


and  X, (c)  =  c,a   +  c^h    ,  X  (c)  =  c^a^  +  C2b2 
rise  to  the  linear  system 

^1^1  ^  ^2^1  -  '^'1  =  ° 

^1^2  ^  ^2^2  -  ^2  =  °  • 


Formulation  I  gives 


(2.7) 


On  the  other  hand,  using  Formulation  II  we  obtain  two  quadratic  equation; 

0 


(c^a^  +  C2b^  -  X^)(c^a2  +  C2b2  -  X^) 
(^1^1  ^  ^2^  -  ^2)^^1^2  -^  =2^2  -  ^2^  =  0 


(2.8) 


If  we  applied  Newton's  method  we  would  find  the  solution  of  (2.7)  in 


one  step,  but  it  would  take  several  steps  to  solve  (2.8)  accurately. 

The  behaviour  observed  in  this  example  occurs  in  general. 

Formulation  II  complicates  the  problem  and  frequently  gives  rise  to 

badly  conditioned  systems.   The  problem  is  intrinsically  simple  along 

some  directions.   Let 

n 
A(a)  =  A(ac)  =   J   ac.A.  =  aA(c)  . 
1=1    ^  ^ 

The  eigenvalues  of  A (a)  are  linear  functions  of  a   and  the  eigenvectors 
are  constant.   Therefore  along  any  1-dimensional  subspace  the  problem 
is  linear.   Formulation  I  preserves  this  property,  which  helps  to 
accelerate  the  convergence  of  the  numerical  methods.   On  the  other  hand 
it  is  easy  to  see  that  using  Formulation  II,  g  is  not  linear  along  any 
1-dimensional  subspace,  but  its  components  are  in  general  polynomials 
of  order  n. 

The  first  mehtod  we  will  consider  consists  of  applying  Newton's 
method  to  solve  the  nonlinear  system  given  by  Formulation  I.   If  all 
the  eigenvalues  of  A(c)  are  distinct  then  they  are  dif f erentiable 
(see  Theorem  1  below)  and 

3X. (c)         „ 
3~     =  q^Cc)'  A.q^(c)  ,  (2.9) 

where  q. (c)  is  an  eigenvector  of  A(c)  corresponding  to  X.  and 
tq.(c)l2  =  1-   As  A(c)  is  symmetric,  ^^q^  ^  1  is  an  orthonormal  set.   We 
can  now  apply  Newton's  method  to  solve  (2.5): 

Method  I 

Choose  a  value  for  c  .   For  k  =  0,1,2,... 

1.  Construct   A(c^)  -    V   cj^  A. 

i4i  ^  ^ 

2.  Find  the  eigenvalues  and  eigenvectors  of  A(c): 

q'^  A(c)Q  =  A 

where  Q  =  [q^  (c  ^)  ,  .  .  .  ,qj^  (c'^)  ]  and  A  =  diag  (X^  (c^)  )  . 

Let  f .  (c^)  =  X  .   (c^)  -  A*   i  =  1,  . . . ,n.   Stop  if 

k         "^ 

I f (c  )!_  is  sufficiently  small. 

3.  Form  J{c^)    by 


J^j  =  qiCc^)"^  A^q^(c'^)  .  (2.10) 


k+1 
4.   Compute  c     by  solving 


J(c'^)  (c^^^  -c^)  =  -f(c^  .  (2.11) 


In  the  case  of  the  additive  inverse  eigenvalue  problem  this  method 
reduces  to  that  of  Downing  and  Householder  C1956).   Even  though  it  is 


stated  in  a  different  form,  it  can  be  seen  that  the  method  used  for 
shell  model  computations  in  nuclear  spectroscopy  (see  Brussard  and 
Glaudemans  (1977))  is  also  equivalent  to  this  method.   Kublanovska ja 
(1970)  gives  conditions  that  guarantee  the  convergence  of  Method  I, 
using  a  line  search  at  step  4. 

Now  let  us  consider  applying  Newton's  method  to  Formulation  II 
Note  that  the  i-th  equation  in  (2.6)  can  be  written  as 

n  * 

g,  (c)  =   n  (X,  (c)  -  A,)  . 


k  =  l 


Therefore,    using     (2.9) 


n      3X,  (c) 

k^i    ^^j 


n     (X, (c)  -  x.) 

Jl?^k        ^  ^ 


I       4    A       q^^       n        (X^-X*) 


k=l 


£7^k 


(2.12) 


Let   G(c)    be    defined    by   G .  .    =    (— ).       Then     (2.12)    can    be    written    as 


G(cj 


n   (X.-x J 


n   (A  -X  ) 


n 


JlT^n 


(A„-X*)  ^ 


^\-\^ 


<{[    A^q^ 


q      A^  q 
^n      l^n 


^1    \   ^] 


q      A      q 
^n      n    ^i 


11  n      1 


=    diag(g.)    diag(j^)W  J 


(2.13) 


where  f^  is  the  i-th  component  of  (2.5),  J  is  given  by  (2.10)  and  the 
matrix  W  is  defined  by 


n   1 


Then  the  Newton  step  for  Formulation  II  is: 


c     =  c   -  G  (c  )    g  (c  ) 


c'^  -  J  -^(c'^)w".^(c^-)  diag(f^(c^))diag(   ^  ^  )g  (c^) 

g^(c  ) 


=  c^  -  J  -"-(c^yw  •'■(c^)f  (c''^)  „ 


(2.14) 


Method  II  differs  then  only  from  Method  I  in  step  4,  where  J(c  )  should 

k     k 
be  replaced  by  W(c  )J(c  ).   This  method  has  been  used  by  Biegler-Konig 

(1981)  and  generalizes  an  algorithm  of  Lancaster  (1964-a) . 

Note  that  if  the  O^*)"   are  all  distinct  then  W -*- I  as  II  f  (c'^)  II  2  ^  0 . 

Asymptotically  (2.11)  and  (2.14)  coincide  in  this  case,  and  both  methods 

are  usually  quadratically  convergent.   Numerical  tests  show  that  Method 

I  converges  faster  (see  the  example  1  in  Section  5)  and  that  Method  II 

frequently  suffers  from  ill-conditioning.   This  can  be  explained  by 

observing  that  W  becomes  nearly  singular  in  various  situations, 

particularly  far  from  the  solution. 


approach  the  solution.   This  destroys  quadratic  convergence.   3y  contra; 
we  show  in  Section  4  below  that  Method  I  is  usually  quadratically 
convergent,  even  in  the  multiple  eigenvalue  case. 

When  the  matrices  {A.}  are  not  sparse  the  cost  of  forming  J  (see 
(2.10))  can  be  very  high.   One  could  avoid  this  problem  by  using 
Broyden's  method  (see  Broyden  (1965)).   However,  numerical  experiments 
indicate  that  Broyden's  method  generally  needs  a  large  number  of 
iterations  to  approach  the  solution. 


3.   Quadratic  conver-^ence  of  Method  I  when  A(c  )  has  distinct 
eigenvalues . 

v;hen  {  ■'  •  }  1  are  distinct  Method  I  has  local  quadratic  convergence 
under  mild  contitions.   This  follows  from  standard  properties  of 
Newton's  method,  using  the  well-known  fact  that  the  eigenvalues  are 
dif ferentiable  when  they  are  distinct. 

Theorem  1 .   Assume  that  (1.1)  has  a  solution  c  ,  i.e.  f(c  )  =  0,  and 
that  A(c*)  has  distinct  eigenvalues.   Then 

(a)  f  is  (Frechet)  dif ferentiable  in  a  neighborhood  of  c  . 

(b)  If  in  addition  the  Jacobian,  which  we  will  henceforth 
denote  by  Df ,  is  nonsingular  at  c*,  then  there  exists  a 
neighborhood  of  the  solution,  N(c*),  such  that  for  all 
c°eN(c*)  the  iterates  of  Method  I  converge  to  c*  and 

-  .       I  c  — c    1 

lim  sup   -5 j— ^   <  " 

k  ^co    Ic  -c  Ot 

In  other  words.  Method  I  has,  locally,  at  least  Q-quadratic  convergence 
in  the  sense  of  Ortega  and  Rheinboldt  (1970) . 

Proof:   The  fact  that  the  eigenvalues  of  A(c  )  are  distinct  implies 
that  the  eigenvectors  are  continuous  at  c   (see  Ortega  (1972),  p.  54). 
The  partial  derivatives  of  f  are  given  by  (2.9)  (see  Ortega  (1972)  or 
Kato  (1966),  p. 81),  and  are  therefore  continuous  at  c*.   It  follows 
that  f  is  dif ferentiable  at  c*.   It  will  be  shown  later  that  the 
partial  second  derivatives  of  f  are  given  by  (4.4),  and  again,  these 
are  continuous,  so  f  is  twice  dif ferentiable  at  c  .   The  quadratic 
convergence  follows  by  the  well-known  properties  of  Newton's  method 
(see  Ortega  and  Rheinboldt  (1970),  p. 312).     D 


4.   Quadratic  Convergence  of  Method  I  when  A(c  )  has  multiple 
eicenvalues . 


Suppose  that  A(c  )  has  multiple  eigenvalues.   Then  the  eigen- 
vectors {q. (c*) }  are  not  unique,  and  they  cannot  in  general  be  defined 
to  be  continuous  functions  of  c  at  c*.»  Furthermore,  the  eigenvalues 
and  hence  f (c)  are  not  in  general  dif ferentiable  at  c  .   These  remarks 
may  be  verified  by  considering  the  example 


"1    0 
0   -1 


0  1' 

1  0 


Note  that  if  c  and  c  are  such  that  A(c)  and  A(c)  have  distinct  eigen- 
values and  Ic  -  c  1,  3c  -  c*I  and  hence  tic  -  c!l   are  arbitrarily  small, 
Iq.  (c)  -q.  (c)  !  can  be  large.   It  is  thus  natural  to  expect  that  Method 
I  has,  at  best,  slow  convergence  to  c  ,  since  standard  proofs  of 
quadratic  convergence  require  f  to  be  dif f erentiable  and  Df  to  be 
Lipschitz  continuous  at  the  solution.   It  is  however  a  remarkable  fact 

that  Method  I  is  usually  quadratically  convergent  to  a  solution  c   even 

* 
when  A(c  )  has  multiple  eigenvalues.   The  reason  is  essentially  that 

the  eigenvalues  can  be  considered  as  twice  dif ferentiable  functions  of 

a  single  variable  along  any  line  passing  through  the  solution  c* . 

Consequently  every  Newton  step  produces  an  excellent  estimate  of  the 

solution,  and  the  iterates  {c  }  converge   quadratically  although 

{Df(c  )}  does  not  converge. 

To  prove  this  we  first  need  the  following  result  found  in 

Rellich  (1953). 

Lemma  1  (Rellich).   Let  B(a)  be  a  real  symmetric  n^n  matrix  function 

of  a  single  variable  a.   Assume  that  B(a)  is  analytic'  in  a  neighborhood 

of  a=  0  ,  i.e.  it  can  be  expanded  in  a  power  series  convergent  in  that 

neighborhood: 

2 
B(a)  =  Bq  +  aB-  +  a  B-  +  .  .  .  . 

Let  o   be  an  eigenvalue  of  B(0)  with  multiplicity  m.   Then  there  exist 
scalar  functions  a-  (a)  and  vector  functions  v^  (a)  ,  i  =  l,...,m.,  all 
analytic  in  a  neighborhood  of  a=  0  ,  such  that  a. (a)  and  v. (a)  are  an 
eigenvalue-eigenvector  pair  of  B(a),  i.e. 

i  =  1 ,  . . . ,m  . 
Furthermore  (v.  (a)}  are  orthonorm.al  and  O'.(O)  =  o,    i  =  l,.,,,m  , 


i  ^   ^  I  we  can  set  a, 
i_u  -  a_j  i 

j^  ^'_(,a)  we  would  need  to  asLxri^;  '-,v-/  =  -,-,/  '-2'^-^'    "  ! --^  i' 


For 
a  and   a  (a)  =  -  a.   If  we 
(a)  we  would  need  to  define 
but  these  are  net  analytic.   However  it  is  clear  that:  in  a  one-sided 
neighborhood  of  a=0,  say   0  £  a  <_  a ,  we  can  assume  that  ^ li^) ^- •  •  ^'^-.i'^)  - 
Lemma  1  shows  that  {a. (a)}  ,  a  representation  of  the  eigenvalues 
of  B(a),  have  derivatives  of  all  orders  in  a  neighborhood  of  a=  0. 
Lemma  2,  which  is  a  modification  of  a  result  of  Lancaster  (1964-b) , 
explicitly  shows  the  form  of  the  first  and  second  derivatives. 


By  a  vector  or  matrix  function  analytic  in  a  we  mean  one  whose 
components  are  analytic  in  a,  i.e.  are  convergent  power  series  in  a 


(applying  this  lerrjri.a  once  for  each  distinct  eigenvalue  of  B(0)) 


I(:i)  =  diag(-.  (a)  )  andV(ci)  =  [v,  (a)  ,  .  .  .  ,  v^  (d)  ]  .   Note  that.  V  (a)  "^V  {a)=  I 


and 


Define 


Z(a)  =  V(a)'^  B(a)V(a)  .  (4.1) 


R(a)  =  V(a)'^  B'(a)V(a)  (4.2) 


with  components  {r..(a)}.   Then  in  a  neighborhood  of  a=0, 

a[{a)    =    r^^(a)  =  v^  (a)  B'Ca)v^(a)  (4.3) 

and  n        ^2 

aV  (a)  =  V.  (a)^  B"(ct)v.  (a)  +  2      I  „    /^^  „    ,    ^ 

(4.4) 


Proof :   Differentiating  (4.1)  and  omitting  the  argument  ct ,  we  obtain 

I'  =  R  +  (V*^)  'bV  +  v'^BV  =  R  +  SZ  -  ZS  (4.5) 

where 

T  I        T 
S  =  (V^)  V  =  -V  V' 

with  components  ts-  ■}.   Equating  the  diagonal  elements  of  (4.5)  then 
gives  (4.3).   Equating  the  off-diagonal  elements  of  (4.5)  we  obtain 

0  =  r^  .  (a)  +  s^.  (a)  (a.  (a)-G^(a)  )  . 

Thus  we  have 

r.^(a)  =0,  i   ^   j,   if  a^ia)    =    ^.{a)  (4.6) 

s..ia)    =    ,.(J2,.(,)  .  i  ,=^  J,   if  a,(a)  ^  a-U)     (4.7) 
Differentiating  (4.2)  we  have 

T  T 

Thus 

o\[a)    =    r:.(a)  =  v.  (a)^  B"(a)v.(a)  +   J^  ^ik^ki"  ^ik^ki  • 


By  B'  (:i)  we  mean  the  matri.N  whose  components  are  obtained  by 
differentiating  the  components  of  B  with  respect  to  a. 


Using  (4.6),  (4.7)  and  the  syruTietry  of  R  we  obtain  (4.4). 


Note  that  if  B(0)  has  multiple  eigenvalues,  say  a.(0)  =  a.(0)  for 
so.Tie  i  5^  j  ,  but  o.(a)  ^   a  .  (a)  for  o.  ^  ^ ,    we  must  have  that  r..  -*  0  as 
a-'O,  because  r.-(a)  (and  a.  (a))  are  continuous  at  a  =  0  by  Lemma  1. 


We  are  now  ready  to  prove  the  main  result. 


Theorem  2.   Assume  that  (1.1)  has  a  solution  c  ,  i.e.  f(c  )  =  0,  and 
suppose  that  A(c*)  has  multiple  eigenvalues.   Let   c    be  the  iterates 
of  Method  I.   Assume  that  for  all  k,  A(c  ;  has  distinct  eigenvalues  so 
that  Df (c  )  exists,  and  assume  that  Df (c  )  is  nonsingular  so  that 
Method  I  is  well  defined.   For  every  vector  u  on  the  unit  sphere,  i.e. 
'u*   =  1'  define 


and  let  (a  )  •  (a) ,  (v  )  .  (a) ,  i  =  1, . . . ,n  be  its  eigenvalues  and  eigen- 
vectors, all  analytic  in  a  neighborhood  of  a  =  0  ,  as  given  by  Lemma  1 
Define  J  (ct)  by 

(J^).j(a)  =  (v^)^(a)^  A.  (v^)^(a)  .  (4.8) 

Now  assume  that,  for  all  u  with  IIUJ2  =  1/  ^ ^^^^    ^^  nonsingular  and 
that 

sup   DJ  (0)~-'-n  <  0=  .  (4.9) 

|UB2=  1    "" 

Then  there  exists  a  neighborhood  of  c  ,  N(c  ),  such  that  for  all 
c  e  N(c*)  the  iterates  of  Method  I  converge  to  c   and 


lim  sup    ^.   ^,   -•  <  "   ,  (4.10) 


i.e.  Method  I  is  at  least  Q-quadratically  convergent. 

Proof:  By  the  remarks  following  Lemma  1  we  can  assume  that  for  all  u 
with  'ul  =  1,  (a  ),  (a)  £..._<  (a  )  (a)°  in  a  one-sided  neighborhood  of 
a  =  0,  say  0  £a  <_ct,  and  hence  that  for  such  a  , 

(a^)^  (a)  =  X^(c*  +  au  )  . 


f  (a)  =  f  (c  +  au) 


A^  (c*  +  au)  -  A* 


A  (c  +  au)  -  A 


(Vn(")-^n 


Now  let  c  be  an  iterate  c   generated  by  Method  I;  we  will  omit  the 
superscript  for  convenience.   Define 


•  c  -  c' 


u  =  g  (c  -  c  ) 


Note  that  SUI2  =  1/  f(c)  =  f  (3)  and  that  Df(c),  which  exists  by 

ons 

c  =  c  -  Df  (c)"-"-  f  (c)  . 
Using  the  fact  that  f(c*)  =  0 ,  we  have 


!c  -  c  t  <  ilDf  (c) 


-1, 


Df (c)  (c  -  c  )  -  f  (c)  +  f (c  ) « 


and  hence 


-1, 


Now  by  Lemma  2  and  the  definition  (4.8) , 


3  =  1 


From  (4.11)  we  therefore  get 


Ic-c  3l!J^(S)  "-3  B3f^(5)  -  f  ^  (S)  +  f  ^  (0  )  il 


The  key  point  of  the  argument  now  appears:  although  f (c)  is  not 

dif ferentiable  at  c  =  c*,  f  (a)  is  dif f erentiable  at  a  =  0  and  we  can 

apply  the  mean  value  theorem  (Ortega  and  Rheinholdt  (1970) ,  p.  70)  to 

give 

l5-c*l  <  IJ  (B)"-^l   sup   If'  (B)  -  f '  (y6)"  6  . 
"        0<Y<1 


Since  (f^j)j_(0)  =  (o  )  .   is  defined  and  given  by  (4.4)  in  Lemma  2,  we 
have,  if  c  is  close  enough  to  c   , 

IS  -c*J  <  23j^(B)"^J  Jf^(0)llB^  (4.12) 

It  is  clear  from  (4.4)  that 


sup    If  (0) 

luJ2=l   "" 


Finally  using  (4.9),  (4.12)  and  the  continuity  of  J  (a)  we  obtain  (4.10) 
This  completes  the  proof  of  the  theorem.     D 

Although  we  have  given  Theorem  2  in  terms  of  the  inverse  eigen- 
value problem,  it  is  clear  that  it  could  be  stated  more  generally.   The 
quadratic  convergence  property  for  Newton's  method  thus  holds  for  any 
function  which  is  not  dif f erentiable  at  the  solution,  but  which  does 
have  the  appropriate  Lipschitz  continuity  of  the  derivative  along  rays 
with  end  points  at  the  solution.   A  different  approach  to  proving 
Theorem  1  would  be  to  generalize  a  theorem  of  Stepleman  (1969) .   His 
result  concerns  the  convergence  of  Newton's  method  for  functions  v;hich 
are  not  dif f erentiable  at  the  solution  but  whose  domain  may  be  divided 
into  a  finite  number  of  regions  on  each  of  which  it  is  dif f erentiable . 
For  examiple,  both  his  result  and  ours  show  that  Newton's  method  has 
quadratic  convergence  when  applied  to  the  function  <ti  (a)  =  ]  sina  |  ,  in 
a  neighborhood  of  a  =  0. 

5.   Numerical  Examples 

v;e  now  present  two  examples  to  show  the  numerical  behaviour  of 
Method  I.   The  first  example  has  distinct  eigenvalues  and  the  second 
has  multiple  eigenvalues  at  the  solution.   We  used  a  VAX  11/780'  at 
the  Courant  Mathematics  and  Computing  Laboratory,  with  approximately 
16  decir.al  digits  of  accuracy  (double  precision  arith^T.etic )  . 

Exam.pl e  1.        (Distinct  Eigenvalues)   Consider  problem  (1.1)  with  n  =  4, 
c*  =  [1,1,1,1]"^   , 

A*  =  .7291977090804996  E+01 

X"  =  .1294937131581000  E+02      A*  =  .4  31651174  5658000  E+02 
3  .  4 


5  4  11 
4  5  11 
114  2 
112       4 


7  6       5 
10       8       7 

8  10       9 
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A   =  I  and  A   =  diag(i).   The  starting  point  was  c   =  [  .  85  ,  .  90 , 1 . 5 , 1 . 3] 
The  results  obtained  using  Method  I  and  Method  II  are  as  follows: 


Iteration 

If! 

2,  Method  I 

Df 

D-/  Method  II 

1 

.337  E+01 

.337  E+01 

2 

.216  E-01 

.412  E+00 

3 

.204  E-04 

.192  E-01 

4 

.949  E-11 

.615  E-04 

5 



.708  E-09 

Both  methods  exhibit  quadratic  convergence  but  Method  I  is  faster.  We 
observed  in  almost  all  our.  tests  that  Method  I  took  less  iterations  and 
had  a  faster  speed  of  convergence  than  Method  II. 

Example  2  (Multiple  eigenvalues).   Consider  problem  (1.1)  with  n  =  6, 

.  *     *     * 

c*  =  [1, . . .,1] ,  X. 


are  defined  as  follows: 
A. 


1,  .  .  .  ,6  , 


k   ^k^k  ■  ^k^k 
where   e,   is  the  k-th  column  of  the  identity  matrix  and 


r  6 1  rs 

3.5  2 

0  2 

0    ,  b^  =   0 

0  0 

_  0  J  |_0 


-4  -4 

0  -3 

2   ,  b^=   5 

3.5  6 

0  J  6 


The  starting  point  was  c   =  I. 9, .9, .9, 1,1, 1.1, 1.1].   Method  I  exhibits 
quadratic  convergence: 


1  .146  E+01 

2  .624  E+00 

3  .492  E+00 

4  .207  E+00 

5  .498  E+00 

6  .984  E-01 

7  .384  E+00 

8  .420  E-01 

9  .250  E-03 

10  .278  E-06 

11  .794  E-13 
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